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We report on experiments that demonstrate dynamical instability in a Bose-Einstein condensate 
at the band-edge of a one-dimensional optical lattice. The instability manifests as rapid depletion 
of the condensate and conversion to a thermal cloud. We consider the collisional processes that can 
occur in such a system, and perform numerical modeling of the experiments using both a mean-field 
and beyond mean-field approach. We compare our numerical results to the experimental data, and 
find that the Gross-Pitaevskii equation is not able to describe this experiment. Our beyond mean- 
field approach, known as the truncated Wigner method, allows us to make quantitative predictions 
for the processes of parametric growth and thermalization that are observed in the laboratory, and 
we find good agreement with the experimental results. 
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I. INTRODUCTION 

Since the realization of dilute-gas Bose-Einstein con- 
densates (BECs) [ij [21 [3] there have been many experi- 
ments on BECs in one-dimensional (ID) optical lattices 
Il[31[7llil[10l[IIl[Il[T3[H[l3[Tl^ 

Neutral atoms in periodic potentials demonstrate a wide 
range of phenomena that are familiar from other physical 
systems. Dilute gas BECs are attractive systems in which 
to study these effects due to the flexibility and control 
of the experimental apparatus, combined with the exis- 
tence of microscopic theories that are computationally 
tractable with various degrees of approximation. 

At the simplest level, atoms in an optical lattice are 
described by Bloch waves familiar from the description 
of electrons in periodic condensed matter systems. At 
sufficiently large densities the interatomic interactions 
result in a rich range of nonlinear phenomena. No- 
tably, the quantum phase transition from a superfluid 
to a Mott-insulator jH [21] and generation of bright gap 
solitons [12 have been observed. When the condensate 
is moving relative to the lattice, the nonlinearity can 
induce a dynamical (or modulational) instability that 
leads to an observed thermalization of the condensate 

[i3[ia[i5i[ia[i7i[i8i[ia[2o]. 

There has been a large amount of theoretical work de- 
voted to instabilities in optical lattices [T9l l2Ql l22l |23[ [241 
[251 IMl lia lia [21 |30l |3ll 132 [3^ 

Two types of instabilities are known to occur: dynam- 
ical and energetic. Energetic (or Landau) instabilities 
exist when the superfluid or condensate is not at a lo- 
cal minimum of energy, and can deplete the condensate 
via dissipative processes such as interactions with a ther- 
mal cloud [32] [33]. Dynamical instabilities manifest as 
a exponential growth of certain modes and are related 



to the process of parametric amplification. A dynami- 
cal instability will rapidly deplete even a pure, isolated 
condensate. This process will often be chaotic and cause 
period doubling ^28| (43j, turbulence and vortex forma- 
tion [27], and a loss of coherence of the condensate. A 
dynamically unstable BEC will eventually thermalize to 
a higher temperature. 

On the other hand, the parametric amplification that 
causes the dynamical instability can be used to gener- 
ate interesting states of the condensate. Hilligs0e and 
M0lmer [44l [45] have shown that momentum and en- 
ergy can be conserved in non-trivial collisions between 
atoms in a ID system with a periodic potential. This 
is in contrast to free space, where such collisions are 
energetically forbidden. This process is analogous to 
phase-matching in optical parametric amplification, and 
has been observed with a Bose-Einstein condensate both 
spontaneously and with seeding by Campbell et al. ^46j . 
Related parametric amplification processes have been ob- 
served in systems with modulated lattices [43 . Para- 
metric generation creates sub-Poissonian number corre- 
lations and quadrature entanglement between different 
modes in optical experiments, and is expected to do the 
same for atomic systems [47 . 

To quantitatively simulate dynamical instabilities in 
a BEC it is necessary to use beyond-mean-field method 
for quantum dynamics. In this paper we will be mak- 
ing use of the truncated Wigner method [48] [49] [50] . 
This method has successfully been previously applied to 
condensates in optical lattices by a number of authors. 
Firstly, Polkovnikov and Wang studied the unstable dy- 
namics of an condensate offset in a combined optical lat- 
tice plus harmonic potential [34 , which lead to damp- 
ened Bloch oscillations. Isella and Ruostekoski investi- 
gated the dynamics of a trapped condensate as an optical 
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lattice was switched on and found number squeezing in 
experiments is limited because adiabatic ramping takes 
a very long time [35l [36] . In another work the same au- 
thors found that quantum fluctuations cause dissipation 
via the dynamical instability in the system [37 . Katz et 
al investigated scattering due to the dynamical insta- 
bility at the band-edge of an optical lattice using a 2D 
model [20 . They observe that the structure of the scat- 
tering halo can deviate strongly from the s-wave halo in 
free space. 

In this paper we report on experiments where we ob- 
serve rapid thermalization of a Bose-Einstein condensate 
at the band-edge of a one-dimensional optical lattice. We 
then model our system using a mean-field approach, but 
demonstrate that we are unable to accurately describe 
the dynamics of our system with the Gross-Pitaevskii 
equation. We then employ the truncated Wigner method 
to simulate beyond mean-field quantum dynamics of the 
system starting from a coherent condensate at zero tem- 
perature. In our simulations, we observe dynamical in- 
stabilities that lead to heating in qualitative agreement 
with the experimental results. We attempt quantitative 
modeling of the experiments and compare our numerical 
results to the experimental data. 

This paper is organized as follows. In Sec. |ll| we de- 
scribe our experiments where a BEG is loaded into a ID 
optical lattice and our observations of the resulting dy- 
namics. Section III models the experiments using the 
mean-field Gross-Pitaevskii equation, and describes the 
physics that can be captured by this approach. In Sec.|IV| 



we describe the nature of the dynamical instability in var- 
ious parameter regimes including those relevant to our 
experiments. We introduce our quantum model of the 
system in Sec. \V\ a nd discuss its advantages and limi- 
tations. Sectior uVTl describes our numerical results and 
studies the instabilities present for systems of different 



dimensionality. In Sec. VII we discuss important exper- 
imental considerations that have not been implemented 
in our simulations and use our results to develop a quan- 
titative estimate of the heating rates that are observed 
in the experiments, before concluding in Sec. |VIII[ 



II. EXPERIMENTAL RESULTS 

A. Experimental Set-Up 

Our experimental set-up has previously been described 
by Mellish et al. pjj. Briefiy, we use a double magnetic- 
optical trap set-up to collect and pre-cool ^^Rb atoms. 
The sample is then loaded into a QUIG magnetic trap, 
from which evaporative cooling results in a Bose-Einstein 
condensate of approximately 1 x 10^ atoms in the \F= 
2, mF=^) hyperfine ground state with a condensate frac- 
tion of approximately 80%. The peak density in the mag- 
netic trap is ~ 1.7 X 10^^ m~^ with an axial trap angular 
frequency of cj^^ ~ 27r x 14.6 Hz and a radial angular 
frequency of o;^ = a;^ ^ 27r x 179 Hz. 



In our experiments a weak moving optical lattice is 
suddenly applied to the condensate at an angle of 63 ± 3° 
to the weak axis of the trap. The lattice potential is gen- 
erated from a single laser of wavelength 780 nm detuned 
4.49 GHz from atomic resonance to minimize sponta- 
neous emission. In the lab frame, the two lasers of wave- 
number Ik/, I = kL ~ 27r/(780 nm) and angular frequency 
difference 5 produce a moving sinusoidal potential 



COS (2ki, • V — 5t) ^ 



(1) 



where Vq depends on the intensity of the lasers and the 
interaction with the atoms. For ease of analysis we define: 



' 2m 



(2) 



where m is the mass of the atoms, and s characterizes the 
strength of the lattice. The frequency difference between 
the beams is chosen so that relative to the lattice the 
atoms are moving with a quasi-momentum correspond- 
ing to the Brillioun zone edge (the Bragg condition). In 
the frame of the lattice, the atoms are moving with mo- 
mentum —h\^L- 



B. Bragg Scattering 

We can gain a basic understanding of the system in the 
dilute, collisionless limit. The single-particle eigenfunc- 
tions of the Hamiltonian without the magnetic trap are 
the spatially periodic Bloch states that are superpositions 
of momentum components differing by 2hkL- Diagonal- 
izing this in the limit 5 <C 1 one finds the Bloch states 
at the edges of the two lowest bands are superpositions 
of the two momentum states (in the frame of the 

lattice) . Denoting the Bloch state with quasi-momentum 
hq in the nth band by n), we have the result: 

\kL,l){t) = l + fc^)-J-fe^) e-»"i« 



V2 



(3) 



and that hu;2 — hcoi = Vo/2. This means that an atom 
suddenly placed in the lattice with definite momentum 
is in a superposition of the above modes, and fur- 
thermore the momentum will Rabi oscillate between 
at a frequency uj2 — uji in a process rather similar to a 
Bloch oscillation. 

By applying the lattice for an appropriate amount of 
time one can transfer any fraction of the condensate from 
one momentum state to the other (Bragg scattering). Al- 
ternatively, by phase shifting the optical lattice at the ap- 
propriate time one can non-adiabatically load a weakly- 
interacting condensate into an eigenstate of the lattice at 
either of the lowest two band edges, as described previ- 
ously in Mellish et al. [TT]. 
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FIG. 1: Typical time-of- flight images of a condensate after it 
has been in an optical lattice at the Brillouin zone boundary 
(Rabi cycling) . Atomic clouds to the right are the atoms that 
have been Bragg scattered (excited), (a) Collisionless regime, 
t = 280 /is (applied to freely expanding atoms 10 ms after 
being released from the trap), after 23 ms total time of flight, 
(b) Nonlinear regime, t = 240 /xs (applied to trapped atoms) 
and 29 ms time of flight. The dashed boxes show typical 
regions of interest used for fitting the data. Also, the clouds 
in (b) are further apart and larger than those in (a) due to an 
additional 16 ms of expansion after the lattice is removed, (c) 
and (d) are x-axis profiles of the absorption images (a) and 
(b) respectively, summed along the vertical axis. 



C. Experiments in the Collisionless Regime 

We present two sets of experiments in this paper. In 
the first, after creating the condensate we turn off the 
magnetic trap and allow the condensate to expand for 
10 ms in order to reduce the effect of interactions. This 
allows for the dissipation of essentially all of the mean- 
field interaction energy. We then suddenly turn on and 
hold the optical lattice potential for a total time fol- 
lowed by a further 13 ms expansion (totalling 23 ms since 
being released from the trap) before imaging. In the 
frame of the lattice, the condensate begins with momen- 
tum —hkL^ and as expected we find that the condensate 
undergoes Rabi oscillations between ^lHUl momentum in 
this frame. 

An absorption image after a lattice evolution time of 
t = 280 /is is provided in Fig. [TJa), where we see roughly 
equal populations in the two momentum modes. To find 
the populations in each mode, we perform a bimodal fit 
to the data at momentum -\-hko and —hko. In the re- 
gion surrounding each momentum state we fit a Gaussian 
distribution for the thermal component and an inverted 
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FIG. 2: Fraction of the population that has momentum > 
in frame of the lattice (i.e. the Bragg scattered atoms) as a 
function of the amount of time the lattice is applied after 10 
ms of free expansion. 



parabola for the condensate. We denote the condensed 
population in the region of interest surrounding the orig- 
inal condensate as A^i and that in the second region of 
interest as N2. Similarly we define the total population 
(condensed and non-condensed) in each region as Pi and 

The oscillation between the momentum states is plot- 
ted in Fig. [2J From the frequency of this Rabi oscillation 
(5.9 ± 0.6 kHz) we estimate the strength of our optical 
lattice, and find that s ~ 3.1 ±0.3. The Rabi oscillation 
is damped due to the finite momentum width of the wave 
function of the initial condensate. An exponentially de- 
caying sinusoid is fitted to the data in Fig. |2] for which 
the damping time is 400 ± 200 /is. We investigate the 
damping due to dephasing theoretically in Sec. |IIIB 



D. Experiments in the Nonlinear Regime 

In the second set of experiments the optical lattice is 
suddenly turned on while the condensate is still confined 
in the magnetic potential. The EEC evolves in the com- 
bined potential for a time t before the lattice, and then 
the magnetic trap, are both rapidly switched off. An ab- 
sorption image is then taken after 29 ms of free expansion 
to observe the resulting momentum distribution. For the 
short range of times t in these experiments there is es- 
sentially no change in the atomic density during the evo- 
lution in the trap, however as can be seen in Fig. [1] the 
momentum distribution undergoes significant evolution. 
The broad momentum distribution indicates that atoms 
have been scattered to a variety of different momentum 
states. The major difference between these experiments 
and those described in Sec. |II C| is that the density is 
much greater while the optical lattice is applied, suggest- 
ing that this effect is caused by nonlinear processes. 

We observe that a fraction of the atoms still undergo 
Rabi oscillation between the ifi/cL momentum modes. 
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FIG. 3: In (a) we plot the fraction of the total population 
that has been Bragg scattered as a function of the time the 
lattice is switched on when the atoms are inside the trap. In 
(b) we plot the fraction of the condensed atoms, found with a 
Gaussian fitting procedure, that have been Bragg scattered. 
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FIG. 4: The non-condensed fraction is plotted as a function 
of the time the lattice is left on inside the trap. 



apparently nonlinear process is the main goal of this pa- 
per. 

III. MEAN FIELD DYNAMICS 

In this section we turn to a common tool of choice 
for modelling BEG experiments near T = 0: the Gross- 
Pitaevskii equation (GPE). The mean- field dynamics de- 
scribed by the GPE have often proven to be in excellent 
agreement with observations in a number of experiments 

A. The Gross-Pitaevskii Model 

The time-dependent Gross-Pitaevskii equation for the 
evolution of the macroscopic wave function '^(r, t) of a 
BEG is 

= -£v^V^(r) + V{r)^{r) + [/o|V^(r)2|V^(r). 

(4) 

The characteristic strength of atomic interactions isUo = 
Anh^as/m where is the s-wave scattering length which 
for ^^Rb we take to be = 100 Bohr radii. The external 
potential V(r, t) in the experiment is 

V{v,t) = "^(u^lx'+u^y+cvy) + VL{r,t), (5) 

which is the sum of the parabolic magnetic potential of 
the trap and the sinusoidal optical potential from Eq. ([T]) . 
The trap frequencies are those measured by the experi- 
ment: (cox.ujy.ujz) = 27r X (14.6,179,179) Hz. 

Our simulations begin with 10^ atoms in the Thomas- 
Fermi ground state, which has a peak density of 
no ^ 1.69 X 10^° m-^. 



In Fig. [3] we plot the oscillating populations of both the 
entire sample, and the condensed modes found by the 
bimodal fitting procedure, finding qualitative agreement 
with previous results in [19 . By fitting an exponen- 
tially decaying sinusoid, we observe a damping time of 
Tp = 450 ± 80 /is for the entire sample. However, we 
find that the condensed modes oscillate for much longer 
with a damping time of tn = 800 ± 100 jas. The fact 
that this is longer than the damping time measured in 
the experiments in the collisionless regime is due to the 
reduced momentum width of the condensate inside the 
trap. In the expansion process the interaction energy of 
the condensate is converted to kinetic energy, resulting 
in a broader momentum distribution. 

The fraction of non-condensed atoms as a function of 
the evolution time in the lattice is plotted in Fig. [4j One 
can see that the non-condensed fraction grows rapidly 
until it saturates at slightly above 50%. This rapid gener- 
ation of the thermal component has been termed 'anoma- 
lous heating' by other authors ^51^. Understanding this 



B. Linear Regime 

The density of the system decreases rapidly after re- 
lease from the trap, and so interactions between the par- 
ticles become negligible in the experiments where the lat- 
tice is applied after 10 ms of free expansion. A simple 
one-body theory can therefore be applied to this system 
[taking Uo\ip^\ = in Eq. Q]. We will now derive an 
analytic approximation to the form of the damped oscil- 
lations. The Rabi frequency of oscillation of an atom is 
dependent on its initial momentum, and the momentum 
width of the condensate gives a spread of Rabi frequencies 
whose sum results in the damping of the full oscillations. 

It is worth noting that the dimensionality of the prob- 
lem can be reduced. Any momentum orthogonal to the 
lattice does not affect the oscillation frequency, and there- 
fore is not a part of our dephasing analysis. Thus, we only 
need to look at the distribution of momenta in the lat- 
tice direction. The distribution is a result of the initial 
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FIG. 5: Comparison of our analytic estimate for the damped 
Rabi oscillations due to dephasing of an non-interacting con- 
densate with a full numerical calculation. In both cases 
we have assumed a Gaussian distribution of initial momenta 
which dephase with each other. Circles represent data from 
numerical simulations, and the line is the analytic result in 
Equation [lO] 



momentum within the trap plus the released interaction 
energy and would be difficult to calculate precisely. We 
introduce a few approximations below. 

To simplify analysis, our initial distribution was taken 
to be Gaussian 



P{k) (X exp 



(6) 



where almost all of our kinetic energy comes from the re- 
leased interaction energy /7int = |/7ono. Due to the high 
aspect ratio of the trap, almost all this energy is dis- 
tributed equally in the two tightly-trapped dimensions. 
Given that the lattice is at angle Q ^ 63° to the weak 
axis of the trap, the conservation of energy requires that 



2m 



sm 
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(7) 



This results in an predicted at ^ 1.23 x 10^ m~^. 

We will now analyze the dynamics that these modes 
will undergo. In a periodic potential, the eigenstates are 
superpositions of momentum states corresponding to the 
same quasi- momentum (i.e. hk ± 2nhkL where n is an 
integer) . We can then find the rate of oscillation between 
two modes of the same quasi- momentum. At the band 
edge, we assume the eigenstates to be superpositions of 
only the momentum modes \^kL). The energy difference 
hcoo between the two states at the band-edge 



2^2 



2m 2' 



(8) 



is accurate for small values of 5, and has a 1% error (com- 
pared to direct diagonalization) for the value of 5 = 3.1 
used in these simulations and experiments. The Rabi 
frequency of the oscillation is uoq. 

Next we will find how the oscillation frequency varies 
with the quasi-momentum. We will assume that the 



eigenstates are superpositions of the modes \/Sk + ki) 
and \/Sk — kL) only. In this basis, we see the energy dif- 
ference between the upper and lower states is (to second 
order in A/c): 



(9) 



We can now use this dispersion relation, and the initial 
conditions of Eq, (|6| to find the distribution as a func- 
tion of CO. Performing the Fourier transform yields the 
population of the negative momentum modes in time 



iV(-) = N^-^ 1 



1 [ujot + arctan(/:/T)/2] 



1/4 



(10) 



where r = sm/16ha1 ^174 /is. Note that the functional 
form of this decay is much slower than exponential, and 
goes as for t > r. To give an idea of the validity 

of the approximation in Eq. ([9|, we also performed a nu- 
merical simulation using the linear Schrodinger equation 
and the same Gaussian initial conditions. Fig. |5] demon- 
strates agreement between the analytical estimate and 
the numerical results. 

To compare with the experiment, we fitted an expo- 
nential decaying sinusoid to the first 300 /is of data in 
Fig. [5j yielding a decay time of 954 /is, compared to a 
measured time of 400 ± 200 /is. 

There are several reasons that might explain the dis- 
crepancy between the theoretical and experimental re- 
sults. Firstly, the initial condition of the condensate 
may be quite different from Gaussian due to nonlinear 
and finite temperature effects. Condensates may begin 
with a small, random center-of-mass momentum, giving 
an apparent dephasing effect shot-to-shot. Alternatively, 
the condensate may have larger momentum width due to 
heating or other effects. 

We have estimated the momentum width necessary for 
an exponentially decaying fit to the derived expression 
to return the same decay time as the experiment. This 
yields <j/c ^ 1.75 x 10^ m~^, which is 41% greater than 
our prediction. Approximately twice the predicted initial 
kinetic energy is required to explain the experimental 
results. 



C. Nonlinear Regime 

We now apply the Gross-Pitaevskii theory to simulate 
the condensate dynamics in the combined magnetic and 
optical potential. We perform a full 3D calculation di- 
rectly matching the experimental parameters described 
earlier. The initial condition for the calculation is gen- 
erated via imaginary time evolution to find the ground 
state containing approximately 1 x 10^ atoms. 

In Fig. |6] we display the momentum distribution af- 
ter the lattice applied to the trapped atoms for (a) 0.5 
ms and (b) 1 ms respectively . The broad, thermal fea- 
tures observed in the experiment are not apparent in the 
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FIG. 6: Population per mode in momentum space is shown 
(a) 0.5 ms and (b) 1 ms after the lattice is applied in a three- 
dimensional GPE simulation. The x-axis is the weak trapping 
axis. The condensate is visible at zero momentum and 2hko, 
but no thermal features are visible. Note the "camera" angle 
is different to the experiment and a logarithmic scale is used. 
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FIG. 7: Total population with positive momentum in the lat- 
tice frame for a full 3D GPE calculation corresponding to 
the in trap experiment showing damped Rabi oscillations. A 
cubic spline (solid line, to guide the eye) and an exponen- 
tial decay (dashed line) were fitted to the simulation results 
(circles). The exponential decay constant is r —1.3 ms. 



simulation. However, small ripples in the condensate mo- 
mentum densities are indications of linear dephasing. We 
observe that the momentum distribution of the conden- 
sate modes broadens slowly in time. The center of mass 
of the two condensates also evolve in time as the conden- 
sates climb the walls of the trapping potential. 

In Fig.[7|we plot the populations with positive velocity 
relative to the moving lattice. An exponentially decaying 
sinusoidal fit results in a decay time of 1.3 ms. We note 
that the exponential decay fits poorly to the numerical 
results. Qualitatively, some features of the shape of the 
decay envelope matches more closely to that derived in 
Eq. (10). In particular, the decay is slower than expo- 



nential for short times and there is a notable phase shift 
(or a change in frequency) reminiscent of the second term 
inside the cosine of Eq. (10). 



In the experiment we were able to separate the evo- 
lution of the condensed atoms from the non-condensed, 
as shown in Fig. |3] (b). The damping of the condensate 
oscillations is at least partly caused by the linear dephas- 
ing. Linear dephasing is present in this simulation, and 
should be comparable to that of the experiment. How- 
ever, we see that the numerical result of 1.3 ms is slower 
than the experimental fit of 800 /is. More favorable agree- 
ment was found in Ref. ^9j where the authors applied a 
similar comparison with their experiment. 

The Gross-Pitaevskii simulation did not generate the 
thermal features found in the experiment. However, 
semiclassical theory can explain the occurrence of ther- 
malization by the existence of dynamical instabilities, as 
shall be discussed below. 



D. Dynamical Instabilities 

The Gross-Pitaevskii equation is a non-linear equation 
that is known to exhibit dynamical instabilities. 

Stationary states of the GPE may, under certain cir- 
cumstances, be unstable points of equilibrium. That is, 
perturbations to the state may grow exponentially in 
time. Eventually, for any realistic initial condition the 
growing perturbations will alter the state significantly. 

The presence of dynamical instabilities can be detected 
by linearizing the GPE about the stationary state, or 
equivalently by solving the Bogoliubov-de Gennes equa- 
tions. Quasi-particle modes with imaginary eigenvalues 
are those that undergo exponential growth and are dy- 
namically unstable. Alternatively, modes can have nega- 
tive eigenvalues indicating energetic (Landau or thermo- 
dynamic) instabilities. 

There has been significant theoretical and experimen- 
tal research into the presence of dynamical instabilities in 
condensates in optical lattices [HI [Ml [T5l [T6l ITT l [22 l [23 l 

ElEailiETlEillglEolEIll^ What 

is certain is that dynamical instabilities cause the ther- 
malization of the condensate when placed at the band 
edge of an optical lattice at high enough densities. 

We will not repeat previous work, but simply point out 
that quantitatively accurate simulations of the regime in 
this experiment have not been performed as far as we are 
aware. This regime involves a condensate at the band- 
edge of a relatively weak ID optical lattice, where trans- 
verse excitations play an important role (as discussed in 
Sec. IV). The reason for this is the lack of spontaneous 



processes in deterministic, semiclassical simulations. Ac- 
curate quantum dynamics in other regimes has been per- 
formed in [Ml Ea EH]. 

The exponential growth of quasi-particle modes is a re- 
sult of Bose enhanced (stimulated) processes. However, 
even with zero initial population in the quasi-particle 
mode there is the chance of spontaneous scattering into 
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that mode. These spontaneous processes are very impor- 
tant when determining the subsequent dynamics of the 
condensate. In fact, in some regimes the spontaneous 
processes dominate over the stimulated ones, as discussed 
in the next section. 



IV. REGIMES OF DYNAMICAL INSTABILITY 
IN ONE-DIMENSIONAL LATTICES 

In this section we discuss the effect of dynamical insta- 
bility on the condensate when one considers the sponta- 
neous processes that may occur. 

The character of the dynamical instability in a lattice 
is dependent on many factors, including the interaction 
strength, the number of particles, the strength of the 
lattice and the dimensionality of the system. The role of 
excitations in the directions transverse to the lattice was 
studied thoroughly in Ref. [30]. They showed that, for 
small enough 5, that there was a significant difference in 
the stability diagram for ID and 3D models. 

In our system, where s 3.1, we would expect trans- 
verse excitations to be important features of the dynam- 
ical instability. Our eventual goal is to quantitatively de- 
scribe the dynamics of the condensate and so eventually 
must perform a full three-dimensional analysis. However, 
we will first discuss the dynamical instability in systems 
of lower dimensionality. 

For quasi-one-dimensional systems, one can expect 
a narrow band of modes for which collisions are both 
resonant (i.e. conserve energy) and conserve momen- 
tum modulo 2hkL (i.e. quasi- momentum in the lattice). 
These are the modes which will undergo rapid exponen- 
tial growth. For the case of a pure condensate, the 
process is initiated by spontaneous collisions between 
condensate particles. Then, as the number of particles 
in the resonant modes grow, Bose-enhancement of the 
scattering process causes parametric amplification. The 
timescale for this growth to become significant will de- 
pend on the initial rate of spontaneous scattering, and so 
should be included in any model which attempts to quan- 
titatively model dynamics. At later times, secondary col- 
lisions between atoms will lead to thermalization of the 
system. 

For higher-dimensional systems there will be a much 
larger volume in phase-space where collisions are close 
to resonance, most of which incorporate an excitation 
transverse to the lattice dimension. In this case the total 
rate of condensate depletion by spontaneous scattering 
is greatly increased. However, as there are many more 
modes the effect of Bose-enhancement is significantly re- 
duced. For some systems the number in each mode may 
not grow much at all before the condensate is entirely 
depleted. In this situation the parametric nature of the 
dynamical instability will be washed out by spontaneous 
processes. Note that secondary collisions would be ex- 
pected in these systems too, which will lead to thermal- 
ization. Thus, in systems of any dimensionality one could 



expect to observe rapid thermalization of the condensate. 

As an example, consider a system in the limit of a 
weak lattice (small s). In this case the dispersion spec- 
trum is not modified significantly from the free particle 
case. However, oscillations in momentum due to the lat- 
tice are expected to occur. At the band-edge one would 
expect Rabi oscillations between modes with momentum 
:^hkL- Assuming that the dispersion spectrum is not sig- 
nificantly different to free space, one would expect reso- 
nant collisions between the modes with momentum ifi/c^ 
into a spherical shell of momentum with modulus hkL 
^2 ^^^3- We can derive the spontaneous collision rate 
using Fermi's second golden rule. 

For a system with a total of ~ 10^ particles in a 
volume V ^ 10~^^ m^ the spontaneous scattering rate 
is: 



dN 



U^N^m\kL\ 



4 X lO^s" 



(11) 



where we have substituted values appropriate to our ex- 
periment. Thus, 10% of the condensate is lost in just 
250 /is by spontaneous collisions. 

Of course, this value will not be accurate for s 3, 
and does not take into account that the atomic popula- 
tions are oscillating between the two momentum states 
or the effect of Bose-enhancement. This does however 
highlight that rapid thermalization does not require para- 
metric growth in dynamically unstable 3D systems, but 
could be simply provided by a sufficiently large density 
of states at resonance. 



V. QUANTUM MODEL 

The second-quantized Hamiltonian for an ultra-cold 
Bose gas is: 



H 



j d^r V^t(r) ^- 



2m 



V2 + y(r,t)U(r) 



(12) 



where the field operator '^(r) annihilates an atom of mass 
m at position r and obeys the commutation relation 

[ViUr),V^(rO] =(^^(r-rO. 

It is not possible to analytically solve for the evolu- 
tion of a Bose gas system directly from the Hamilto- 
nian in Eq. (12). Analytical approaches always involve 



approximations, and often linearize the quantum fiuc- 
tuations about the mean-field as in the Bogoliubov ap- 
proach. Numerical solutions are possible with the exact 
Positive-P method, but quickly become intractable with 
sampling problems in many situations [48 . The approxi- 
mate Gross-Pitaevskii equation (GPE) is straightforward 
to solve, but fails to account for any spontaneous effects 
that are important for the problem we are considering. 

In this paper we implement the truncated Wigner 
phase-space method for the system dynamics. This is 
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an approximate method that goes beyond the GPE by 
including both spontaneous processes and the dynamics 
of non-condensate modes. 



A. The Truncated Wigner Method 

We will now briefly explain the origin and use of the 
truncated Wigner method. This was used in the context 
of squeezing of solitons in optical fibres by Carter and 
Drummond [56 , and was first applied to Bose gases by 
Steel et al [48 . An equation of motion for the Wigner 
function can be obtained from the Hamiltonian Eq. (12) 



by deriving the master equation for the system followed 
by the appropriate transformation of the density oper- 
ator. However, for a system with n modes the Wigner 
function is 2n-dimensional and it is impractical to solve 
such an equation directly. 

However, truncating the third-order derivatives in the 
equation of motion for the Wigner function results in 
a Fokker-Planck equation, which can then be simulated 
using stochastic methods [49^ ^ . The omission of the 
higher order terms has previously been justified when the 
there are more particles than modes in a calculation and 
the simulated times are short [54, 55l The equa- 

tion of motion for the stochastic trajectories is simply 
the usual Gross-Pitaevskii equation, although the initial 
conditions for the field are sampled stochastically from 
the appropriate Wigner distribution. A large number of 
trajectories must be calculated, which are then sampled 
to obtained expectation values of symmetrically-ordered 
quantum field operators. 

As an example, the ensemble average is equal to 
the expectation of the symmetric ordering of the num- 
ber operator, or {iIj^iIj + iljiilj^)/2. Therefore, the average 
number of particles in a mode is (i^^i^) = — 1/2. For 
more details on the truncated Wigner and other phase 
space methods see [48"^ "JS, 50 . 

Because of the nonzero width of the Wigner distribu- 
tion, every vacuum mode begins with uncorrelated Gaus- 
sian noise, normalized to an average of half a particle per 
mode. This provides the seeding needed for the equation 
of motion (the GPE) to allow scattering events into these 
modes. Furthermore, the rate of scattering corresponds 
to that of the expected spontaneous processes. 

Allowing such spontaneous processes is essential to ac- 
curately model systems exhibiting instability. Unstable 
dynamical solutions of the GPE do not correspond accu- 
rately to quantum mechanics. Such states would eventu- 
ally be altered by spontaneous processes and are therefore 
not equilibrium states. The truncated Wigner approach 
unambiguously allows both dynamical instabilities and 
Landau instabilities to manifest. Additionally, a back- 
ground thermal gas can be included in the initial condi- 
tion to allow dissipation via the Landau instability. The 
importance of the two instabilities has been a point of 
discussion in the past [13 . 

The disadvantage of the truncated Wigner method is 



the accuracy of its fundamental approximation is difficult 
to quantify. However, some work has been done which 
demonstrates its validity regimes [55] [58l [59l |60l [61]. 
Some of these works have concluded that the major re- 
quirement for accuracy is that the number of modes 
must be much smaller than the number of real parti- 
cles [58l [59]. If this is not the case, then the "vacuum 
noise" in the initial condition can no longer be consid- 
ered a perturbation. The vacuum fluctuations will begin 
to interact as if they are real particles, and unrealistic 
effects such as negative expectation value for population 
become prolific [59] . 

The method in [60l [61] allows one to check the validity 
of the truncated Wigner approximation at the expense 
of significant extra computation. In this paper, however, 
we instead choose to work in the regime where there are 
more particles than simulated modes. 

Unfortunately, this is not the case for the full 3D exper- 
imental parameters. Including the trap and the optical 
lattice whilst simultaneously avoiding numerical aliasing 
due to the nonlinearity requires of order 10^ modes to 
represent 1 x 10^ atoms. Our calculations in this regime 
clearly display unphysical behavior similar to those found 
in [59 . In particular, we observe large areas of negative 
population at high momentum modes, and a correspond- 
ing clustering of the vacuum noise at small momentum 
which swamps the condensate. The condensate is rapidly 
depleted as it scatters off these unphysical particles. 

To proceed we must make a further simplification by 
neglecting the trapping potential, and instead modeling 
a homogeneous system with a similar density to the ex- 
periment. By doing so we expect to capture the most 
important physics, as the timescale that the thermaliza- 
tion manifests is much shorter than the inverse trapping 
frequencies, so that whilst the optical lattice is on the 
atoms will be essentially stationary and the density con- 
stant. The most important effect of the trap in the exper- 
iment is to maintain the high density so that collisions 
may occur. This model requires fewer modes to simu- 
late the dynamics, and therefore the truncated Wigner 
method is expected to be more accurate. 

In the next section we describe the details and present 
the results from our truncated Wigner simulations. We 
point out in advance that while these calculations will 
model the spontaneous scattering in this system, the ho- 
mogeneous calculations will not capture the linear de- 
phasing of the Rabi oscillations in the experiment as the 
condensate begins in a single momentum state. We will 
consider this further in Sec. I VIII 



VI. TRUNCATED WIGNER SIMULATIONS 

In the experiment the initial inhomogeneous conden- 
sate wave function spreads over about 27 wells of optical 
lattice potential. We use a total population of 10^ atoms 
and choose to use a density of no = 10^^ m~^, roughly 
corresponding the the average density within the trap 
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rather than the peak density. We have implemented a 
rectangular grid with periodic boundary conditions for 
our simulation space, with dimensions x Ly x L^. The 
lattice is parallel to the long axis of the space, with spac- 
ing 390 nm, extending over 32 wells, making = 12.48 



/im. It follows that Ly 



8.95 /im. Note that 



this geometry is slightly different to that of the extended 
condensate at a 63° angle to the lattice. 

We begin with a homogeneous cloud moving at the 
Bragg condition with momentum of hkL in the direction 
of the stationary lattice. The nonlinearity and strength 
of the lattice are the same as used previously. Care was 
taken to ensure numerical accuracy, and a projection op- 
erator that prevents numerical aliasing [62l |63l |64] was 
employed. 

In the following subsections we present our results of 
the truncated Wigner method used to simulate the ex- 
periment using ID, 2D and 3D models. This was done 
for two reasons. The first was to demonstrate the differ- 
ent regimes discussed in Sec.[lVl Secondly, the truncated 
Wigner method will be very reliable in the ID and 2D 
models where the number of modes is much smaller than 
the number of particles. 

To compare the simulations at different densities we 
have ensured that the nonlinearity has scaled correctly. 
In ID, the nonlinearity is Uq/ {LyLz) and in 2D it is 



The one-dimensional results in Sec. VIA clearly show 
the dynamical instability and the parametric gain of the 
resonant modes. The two-dimensional simulations in 



Sec. VI B allow us to clearly visualize the transverse ex- 
citations, as was achieved previously in [20 . We also 
analyze the coherence and statistics of the field in 2D 
where statistical errors are more manageable that the 3D 
situation, as more trajectories can be computed for the 
same amount of CPU time. Finally, in Sec. | VI C] we show 
the results of the complete three-dimensional simulation 
that exhibits strong spontaneous scattering and make a 
comparison to the experimental results. 



A. One Dimension 

We first investigate the dynamical instability in one- 
dimension. Figure [8] shows results for a single trajectory, 
which in some sense can be considered to be analogous to 
what might be observed in a single experiment in the lab 
[55] . In real space, we observe spontaneous breaking of 
the translational symmetry of the lattice brought about 
by the dynamical instability. The high-density features 
correspond to generation of bright-gap solitons. These 
solitons are known to be generated from states near the 
top of the ground band for condensates with repulsive 
interactions [23 . The momentum distribution broad- 
ens indicating spontaneous collisions between condensate 
atoms scattering atoms to other momentum modes. 

The work of Hilligs0e and M0lmer [44[ |45] predicts 
that collisions that conserve energy and momentum (i.e. 
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FIG. 8: (Color online) Top: Density of a ID BEC moving 
in an optical lattice. We see spontaneous symmetry breaking 
of the BEC in in a single trajectory simulation. The local- 
ized regions of high density are bright-gap solitons. Bottom: 
The same results are presented in momentum space, where 
symmetry breaking is characterized by mixing into a variety 
of different modes (note the logarithmic dependence of the 
density scale). 



resonant collisions in the eigenbasis of the lattice) will ex- 
perience parametric growth. If energy conservation is not 
quite satisfied for a particular mode then observe Rabi 
oscillations are expected that will hinder the exponential 
growth of parametric amplification. This is exactly what 
we observe in the simulation as the resonance condition 
is never precisely matched. This effect is demonstrated 
in Fig. |9j where population revival into the states 
can be seen. Also note that only about half the popu- 
lation is depleted from the two condensate modes - this 
corresponds to the half that are in the lowest band near 
the band-edge that generate a bright-gap soliton [23^. 

The experiments described in this paper indicate 
damping of the Rabi oscillations is due to dynamical in- 
stability. Our ID simulations show damping occurring 
after a few milliseconds (see Fig. |9|. Ensemble averages 
of 1000 trajectories have shown the spread of momen- 
tum states and the approach towards thermalization. It 
should be noted that phase coherence is also destroyed by 
the dynamical instability. This is an interesting example 
of decoherence under Hamiltonian evolution. 
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FIG. 9: Top: Sum of the population in the momentum modes 
zLkL of a ID Bose gas plotted versus time. The small, rapid 
oscillations are due to the initial state mixing with higher 
band states in the Bloch basis. The atoms are depleted from 
the condensate mode, until a partial revival begins at ~ 2ms. 
This is due to Rabi oscillations of the not-quite resonant col- 
lisions. Bottom: Total populations of atoms with positive or 
negative momentum versus time. Damping of the Rabi oscil- 
lations is observed in the ID simulations after 1.5ms, but then 
revive at the same time as the condensate as shown in the top 
figure. Both plots have used an ensemble of 1000 trajectories. 



B. Two Dimensions 

1. Dynamical instability 

As described above, the effect of dimensionality on the 
system varies according to the parameters of the exper- 
iment, and has been dealt with analytically in detail in 
[30 . In that treatment, the dynamical instability was 
studied as an exponential growth of Bogoliubov modes. 
A complete set of Bogoliubov modes includes excitations 
in the direction of the lattice (such as would be found 
in a ID treatment), excitations in the directions perpen- 
dicular to the lattice, and excitations which are mixtures 
of both. For some parameter regimes, the perpendicular 
excitations are an important part of the dynamics. 

The same is true in our homogeneous treatment. Col- 
lisions may be resonant into modes with non-zero com- 
ponents of momentum in the directions perpendicular to 
the lattice. This can be observed quite clearly in the 
results of our two dimensional simulations. 



Figure 10 shows a comparison momentum distributions 
for four different initial conditions averaged over 10^ tra- 
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FIG. 10: (Color online) Populations in momentum space of 
a 2D simulation, (a) The condensate began in the lattice in 
the —kL momentum state, t = 1.1 ms. (b) The condensate 
began in a superposition of zL/cl momentum in free-space for 
comparison with (a), t = 1.1 ms. (c) The condensate be- 
gan in the lowest band-edge of the lattice, t = 0.36 ms. (d) 
The condensate began in the edge of the first excited band, 
t = 1.1 ms. Movies of these s imulations are available online 
at: http:/ / www. physics . uq. edu . au/people / f erris / ^ 
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FIG. 11: (Color online) Plots of modes for which the reso- 
nance condition is satisfied. The dotted lines correspond to 
regions close to resonance, having an energy mismatch cor- 
responding to a Rabi oscillation period of 1 ms. The large 
pluses represent the initial state of the BEC. (a) Two atoms 
colliding within the lowest band, (b) Two atoms begin in the 
second band-edge can collide to result in an atom in the low- 
est and second band (inner circles, black online) or two atoms 
in the lowest band (outer rings, red online). 



jectories. In Fig. 10 ^a) the condensate begins in the 
momentum mode; in Fig. 10 'b) we make a comparison 
to free space collisions between momentum modes i/c^. 
Figs. [lOjc) and (d) begin with the atoms loaded into the 
ground and first excited band-edge states. 

When the moving lattice is switched on such that the 
BEC is in the -\-kL momentum state, we see complicated 
dynamics including Rabi oscillations and scattering into 
a variety of modes. A significant percentage of the atoms 
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have developed momentum in the direction perpendicular 
to the lattice, so we obviously can not characterize the 
system as quasi- ID. 

We can understand the complicated dynamics by real- 
izing that the initial condition is a superposition of mo- 
mentum states as given in Eq. (|3|. Each of these states 
has a well-defined energy, and we can analyze which col- 
lisions will conserve energy and quasi- momentum (mod- 
ulo 2kL) in the Bloch basis. The system has been ana- 
lyzed previously in a similar fashion [20l [45] . Katz et al. 
[20] provided comparisons with the results of truncated 
Wigner simulations and an experiment. They found, as 
we do, that the structure of the scattering halo differs 
significantly from the s-wave scattering sphere generated 
by condensate collisions in free space, due to the modified 
dispersion curve. 

Two atoms initially at the edge of the first band can 
collide to produce atoms elsewhere in that Bloch band, 
but with additional kinetic energy in the transverse direc- 
tion. We have calculated the resonance condition using 
the Hartree-Fock mean field method f65^. The values of 



-^coherent (k, t) -^incoherent (k, t) 



momentum where this occurs is shown in Fig. 11 'a) (c/. 
Fig.jlOj^c)). X There are two possible outcomes for atoms 
colliding at the first excited band edge. Either one or 
both atoms can end up in the lowest band. Both of these 
cases are presented in Fig. [iTJb) [cf. Fig. [Toj^d)]. Note 
that the resonance condition is sharper than in Fig 11 'a), 
i.e. there is a smaller area of modes which have a final 
energy close to the resonance condition. This explains 
the results in Fig. [lO] where the lowest band-edge initial 
condition exhibits significantly more scattering at short 
times than the second band-edge case. 



2. Coherent and incoherent components 

We have investigated the phase coherence of the atomic 
cloud as it undergoes thermalization. A Bose-Einstein 
condensate has off-diagonal, long-range order where there 
is a well-defined relative phase across the system. Any 
overall phase of the condensate results from spontaneous 
symmetry breaking, and the value of the global phase is 
undetermined. 

Nevertheless, when employing classical field methods 
one typically begins with a coherent state, ascribing a 
global phase to the condensate at t = such that: 



(13) 



where is the annihilation operator for the condensate 
mode. All observables of the closed system being simu- 
lated will return the same results if one averages over en- 
sembles of random global phase, that is averaging arg(aQ) 
from to 27r. 

The value of the phase of (a^) is only meaningful in 
respect to a second condensate of atoms. An interfer- 
ence experiment can resolve the phase relationship of 
two condensates. If one condensate undergoes heating 
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FIG. 12: (Color online) Coherent (a,c,e,g) and incoherent 
(b,d,f,h) populations in momentum space for the 2D simu- 
lation at times (a,b), 0.56 ms (c,d), 1.1 ms (e,f), and 3 ms 
(g,h). 



and becomes thermalized in an experiment, then that 
cloud of atoms will lose its coherence and not result in 
fringes in an interference experiment. For a thermal sam- 
ple (a) = 0. 

We define the coherent component of the cloud to be 
that with well-defined global phase 

A^coherent(k)= |(a(k))|'. (14) 

The incoherent component is the remainder of the atoms 

^^coherent(k). (15) 



A^incoherent(k) = (a'^(k)a(k)) 



We plot the coherent and incoherent populations in 



Fig. 12 for the initial condition with the condensate in 
momentum state l+Zc^). Coherent dynamics caused by 
the lattice drives the condensate into the modes with 
momentum i/c^;,, iS/c^;,, etc. Only these modes contain 
any significant coherent population. Note the noise in 
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FIG. 13: (Color online) The normalized second order 
momentum-space correlation function g^"^^ (k) is shown for the 
2D simulation at times (a) ms, (b) 0.56 ms, (c) 1.1 ms, and 
(d) 3.0 ms. Results are only shown for modes with average 
population greater than 1/2. 
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FIG. 14: Time dependence of the total population of all 
modes with positive or negative momentum in the lattice di- 
rection (in the frame of the lattice). These exhibit damped 
Rabi oscillations in the 2D simulations similar to the ID re- 
sults in Fig. [9] Note that onset of damping is more rapid in 
2D than ID, and that revivals are not seen in this case. 



condensate modes continue to display coherent statistics. 
However, we see that the all other populated modes have 



/2) 

which supports the claim that a thermal cloud is grow- 



2. These modes are displaying Gaussian statistics 



ing around the condensate. Comparing Figs. 12 and 13 
we see the modes with phase coherence display second- 
order coherent statistics and those without phase coher- 
ence display thermal statistics. Note again that the noise 
in Fig. [13] is statistical. 



Fig. 12 (e,g) is a statistical artefact from averaging over 
only 1000 trajectories, and could be reduced with more 
CPU time. 

The incoherent population is generated by spontaneous 
scattering into the non-condensate modes surrounding 
the condensate. At short times energy conserving colli- 
sion from the two condensate modes dominates. At later 
times further scattering from the newly populated modes 
broadens the distribution and results in thermalization. 



3. Local correlation functions 

The normalized second-order local momentum- space 
correlation function is: 



4' Rabi oscillations 



(at(k)at(k)d(k)a(k)) 
(ank)a(k))^ ' 



(16) 



and can be calculated using the truncated Wigner 
method by transforming into symmetric ordering. This 
correlation function allows us to probe the quantum 
statistics of each momentum mode - in particular, the 
occupation statistics in each mode. For coherent, Poisso- 
nian statistics ^^^^ = 1, and for thermal, Gaussian statis- 
tics ^(2) = 2. 

We plot ^^^^(k, t) in Fig. 13 for the state beginning 
with momentum kL- At all times g^'^\-^kL) ^ 1, and the 



Finally in Fig. 14 we present results the time depen- 
dence of the positive and negative momentum compo- 
nents for the 2D system. These show damped Rabi os- 
cillations with no revivals as observed in in the ID sim- 
ulations. There are many more modes in the 2D system, 
which has two important effects. First, each mode will 
Rabi oscillate at a different frequency, and so one ex- 
pects dephasing to wash out the revival. Second, these 
modes will interact with each other, and the system will 
approach thermal equilibrium faster. We also see that 
these the oscillations undergo faster damping compared 
to the ID results. 



C. Three Dimensions 

1. Numerical accuracy 

The truncated Wigner method for quantum dynamics 
is an approximation that is only known to be accurate in 
certain regimes [55 , 58 , 59 , 6OJ. It is known to be accurate 
under the condition that there are many more particles 
in the simulation than modes. For the experiment we are 
considering this is difficult to achieve, as in three dimen- 
sions the number of modes can easily be greater than one 
million. Thus, in our calculations we took great care to 
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FIG. 15: Condensate occupation of the modes zt/i/ci, versus 
time for a condensate collision in free space. We compare pre- 
dictions from truncated Wigner simulations (solid) , a straight 
line fit to the rate given by Fermi's golden rule (dotted), and 
the solution to the differential equation Eq. ( 11 ) that includes 
condensate depletion (dashed). 



use the minimum possible number of modes without mis- 
representing the physics. To this end we incorporated a 
projector to implement a spherical cut-off in momentum 
space at a radius slightly greater than where the atoms 
can be expected to be found as estimated from lower- 
dimensional simulations. A sufficiently large grid was 
used in real space to eliminate the effects of numerical 
aliasing [62l [63l [64] . This reduces the number of modes 
to ^ 4 X 10^, for the simulation of 10^ atoms. The trun- 
cated Wigner method will not be accurate for long time 
scales under these conditions, but we are confident that 
the dynamics we see in relatively short time scales here 
is accurate. 

We have not implemented the method presented by 
Polkovnikov pOl |6T] to check the validity of our truncated 
Wigner simulation. Instead, we have tested our accuracy 
by comparing our simulation of a simple situation with 
known results. We can compare the spontaneous scatter- 
ing rates from a simulation with atoms with momentum 
with Fermi's second golden rule, Eq. ( 11 ). We have 



performed this simulation for the same parameters, but 
without the lattice. The number of atoms remaining in 
the zb/i/^L momentum modes is plotted as a function of 
time in Fig. 15 At small times the results agree well. 



We can see the effect of Bose-enhancement at later times 
when Fermi's golden rule underestimates the depletion 
rate. 



2. Momentum distribution 

In the experiment, time-of- flight expansion was used to 
image the momentum distributions of the condensates. 



giving 2D column densities. Fig. [16] shows the results of 
the three-dimensional calculation where we have summed 
the populations over the z momentum component. The 
ensemble average is for only 32 trajectories due to the 
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FIG. 16: (Color online) The momentum space column densi- 
ties for the 3D simulation is shown at times (a) ms, (b) 0.56 
ms, (c) 1.1 ms, and (d) 3.0 ms. The scattering halo is clearly 
visible in (c) and (d). 
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FIG. 17: (Color online) The average population per mode is 
shown as a function of axial and radial momenta, for a 3D 
simulation at time 3.0 ms. 



computational demands of the simulation. Note that the 
sharp features seen in the 2D case (Fig. 



10) are some- 



what washed out, resulting in a distribution more closely 
resembling two (compact) condensates surrounded by a 
thermal cloud. However, this is not precisely the case 
- the population is slightly concentrated in ring-shaped 
structures in 3D as seen in Fig. 17 These ring-shape 
structures are similar to those seen in the 2D simulations. 

We remind the reader that the simulations performed 
here are for a homogeneous gas, whereas the condensates 
in the experiment has a finite momentum width. It seems 
reasonable to expect that the sharper features would not 
be as visible in a trapped condensate. In experiments, 
the expansion is done in finite time which can cause sig- 
nificant broadening of the condensate peaks and other 
features. Thus, it is reasonable to say there is good qual- 



itative agreement between Fig. 16 and the experimental 
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FIG. 18: Population of the momentum modes with kx < for 
the 3D simulations. We see damped Rabi oscillations similar 
to Figs.|9]and 14 A decaying exponential fit (dotted line) has 
been made with the numerical results (points), with period 
177 fis and decay time constant 2.00 ms. 
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FIG. 19: The total coherent population of the condensate 
beginning at temperatures corresponding to pure condensate, 
90% condensate, and 80% condensate. All simulations contain 
a total population of 10^ atoms. 



images in Fig.[T]when one takes the effects of the trapping 
potential into account. 



3. Rabi oscillations 

We have plotted the total population with a posi- 
tive component of momentum in the lattice direction in 



Fig. 18 An exponentially decaying oscillation of the form 



Pq cos{u;t)e' 



-t/r 



(17) 



was fitted to the data using a least-squares method. The 
oscillation period was found to be 27r/a; = 177 /is, as 
expected for a lattice with s = 3.1. 

The decay constant was found to be r = 2.00 ms. 
This decay is significantly slower than in the experiment, 
where the decaying fit yields r = 450 /is. There are 
a number of important experimental considerations that 
must be accounted for in addition to this result. We shall 
discuss how to estimate the effects of finite temperature 



and the momentum spread of the condensate in Sec. VII 



resonant collisions. We have performed a 3D simulation 
with the same parameters as described in Sec. |VI C[ but 
with a 10% and 20% thermal fraction surrounding the 
pure condensate centered with momentum -\4ikL- Note 
that the experiments began with approximately a 20% 
thermal fraction. 



In Fig. 19 we plot condensate fractions as a function of 
time. Fermi's golden rule for scattering predicts that the 
scattering rate is proportional to the square of the num- 
ber of condensate atoms (c/. Eq. (pT])). However, such 
a calculation does not take into account the population 
of the modes to which atoms are scattered into. A closer 



VII. EXPERIMENTAL CONSIDERATIONS 



look at the results in Fig. 19 indicate that the fractional 
rate of loss of the coherent component is actually slightly 
greater in the cases where a thermal cloud is present, af- 
ter a time of 200 /is. The Bose-enhancement due to the 
thermal fraction outweighs the reduction in condensate 
fraction. 

The oscillations are also damped more quickly at 
higher temperatures. An oscillating exponential decay 
yields damping times of r ^ 2.00 ms, 1.76 ms, and 1.61 
ms for 0%, 10% and 20% thermal fraction respectively. 
The presence of a thermal fraction can partially account 
for the faster rates seen in the experiment. 



In this section we consider the effects of the physics 
that our homogeneous calculations were not able to cap- 
ture, and estimate their effect on the damping of the Rabi 
oscillations. 



A. Finite temperature 

The simulations we have presented so far have all be- 
gun with initial states sampled from the zero tempera- 
ture Wigner distribution for the condensate. However, if 
the system is in fact at any non-zero temperature then 
any surrounding thermal cloud is able to increase the 
rate of depletion by means of Bose-enhancement of the 



B. Peak-Density Calculation 

Until now our choice of density for the calculations has 
been the approximate average density of the condensate 
in the experiment. Here we repeat the simulations us- 
ing instead the peak density of 1.7 x 10^^ atoms m~^. 
We have no a priori reason to prefer one over the other. 
However, as the scattering will occur quickest in the re- 
gion of highest density, one might expect the dynamics 
in the peak-density region to dominate. 

We performed simulations at zero temperature and 
with a 20% thermal fraction, finding damping rates of 
1.22 ms and 1.00 ms respectively. These numbers are 
closer to the experimental results. In the next section we 
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Description 


Density 


Thermal 


Damping 




(m-3) 


fraction 


time 


Linear Theory 




_ 


954 /is 


Linear Exp. 


_ 


_ 


400 =b 200 /is 


Truncated Wigner 


1.0 X 10^° 


0% 


2.00 ms 


Truncated Wigner 


1.0 X 10^° 


10% 


1.76 ms 


Truncated Wigner 


1.0 X 10^° 


20% 


1.61 ms 


Truncated Wigner 


1.7 X 10^° 


0% 


1.22 ms 


Truncated Wigner 


1.7 X 10^° 


20% 


1.00 ms 


GPE in trap 






1.3 ms 


Combined Theory 


1.7 X 10^° 


20% 


560 /is 


Trapped Exp. 




-20% 


450 ± 80 /is 



TABLE L The rates of damping in the simulations and ex- 
periments. 

combine this with linear dephasing for a full comparison. 

C. Dephasing and Instability 

A summary of the results of all our 3D calculations and 
the experimental results are shown in Table [l[ We can 
now combine the results of the linear dephasing and the 
simulation of the dynamical instability. In the simplest 
model for decay of the Rabi oscillations, we would simply 
add the two exponential decay rates by = T^phasing + 

''^ns\abiiity • dcphasiug fouud in the Gross- 

Pitaevskii model and instability at the peak density of 
1.7 X 10^^ m~^ with a 20% thermal fraction. This results 
in a decay constant of r 560 /is, which is slightly longer 
than in the experiment. Given the uncertainties involved 
in calculating this result we consider that it is in rather 
good agreement. 

VIII. CONCLUSIONS 

In summary, we have described experiments demon- 
strating the interaction-induced instability of a BEG at 



the band-edge of a ID optical lattice. We have quali- 
tatively studied the effects of dynamical instabilities in 
one-, two- and three-dimensional systems, demonstrat- 
ing different regimes of the dynamical instability. We 
have attempted to quantitatively model the dynamics of 
this system by using the truncated Wigner method, with 
good agreement between our experimental and theoreti- 
cal results. 

Our truncated Wigner simulations produced slightly 
slower scattering and damping of Rabi oscillations than 
observed in the experiment and showed that the scat- 
tering rate is larger for non-zero temperatures. How- 
ever, in our simulations we have neglected the trapping 
potential in order find a regime in which the truncated 
Wigner technique was accurate. We have investigated 
some of the consequences of the trap using simple linear 
and mean-field calculations, and found their contribution 
to damping of Rabi oscillations. 

Although the truncated Wigner method has had lim- 
itations in the particular experiments presented here, 
there are many other experiments to which it could be 
applied. In particular, systems of reduced dimensionality 
where the condensate is tightly confined in one or two 
directions (i.e. waveguides) will be well modeled with 
the truncated Wigner method. This will allow quantita- 
tive analysis of non-classical effects in condensates, such 
as spontaneous scattering, number squeezing effects and 
entanglement between atoms and with light. 
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